load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;**************************************
begin

season = getenv("season")
VAR_measure = getenv("VAR_measure")
ilevel = toint(getenv("ilevel"))

season = "JJA"
VAR_measure = "VAR_R95_R50"
ilevel = 4

VAR_measure_regional = VAR_measure+"_regional"


;---------------------

Nmodel = 60
N1 = 35        ;CMIP5
N2 = 25        ;CMIP6


nlat = 180
nlon = 360

Nindex = 4
index_name = (/"R90","R95","R98","R99"/)

Nvar = 3
var_name = (/"pr1d","pr5d","pr30d"/)

Nlevel = 6
warming_levels = (/"0C","0.5degC","1degC","2degC","3degC","4degC"/)
warming_levels_plot = (/"0C","0.5degC","1degC","2degC","3~S~o~N~C","4~S~o~N~C"/)
;level_idx := (/4,5/)

Nobs = 4
dataset_obs = (/"GPCP","GPCC","CPC ","ERA5"/)

NSMILEs = 5
SMILEs_model = (/"CESM1-CAM5","CSIRO-Mk3-6-0","CanESM2","EC-EARTH","GFDL-CM3"/)


;*******************************************
;regional
;*******************************************
lonW := (/0,0,80,115,0,260/)
lonE := (/360,360,150,180,25,300/)
latS := (/45,-70,45,30,50,50/)
latN := (/70,-45,70,45,70,70/)

subregions := (/"NH_extratropics","SH_extratropics","Northern_Asia","WNP_monsoon","Europe","Eastern_NAmerica"/)
subregions_plot := (/"NH mid-latitudes","SH mid-latitudes","Northern Asia","West North Pacific","Europe","Eastern North America"/)
domain := (/"ocean","ocean","land","ocean","land","land"/)
Nregion := dimsizes(lonW)



;--------------read change in frq---------------
diri := "/WORK/zhangwx/EConstraint/data/"
change := new((/Nvar,Nmodel,Nindex,nlat,nlon/),float)

do ivar = 0,Nvar-1

  f1 := addfile(diri+"CMIP5/pr/"+var_name(ivar)+"/"+season+"/RX_pct/multimodel/change_frequency_RX_pct_"+var_name(ivar)+"_multimodel_warming_levels_"+season+".nc","r")
  change(ivar,0:N1-1,:,:,:) = f1->change_frq_RX_pct(:,ilevel,:,:,:)               ;(35model, 4index, lat, lon)

  f2 := addfile(diri+"CMIP6/pr/"+var_name(ivar)+"/"+season+"/RX_pct/multimodel/change_frequency_RX_pct_"+var_name(ivar)+"_multimodel_warming_levels_"+season+".nc","r")
  change(ivar,N1:N1+N2-1,:,:,:) = f2->change_frq_RX_pct(:,ilevel,:,:,:)        ;(25model 4index lat lon)

end do

PR := change/100.+1             ;% change --> probability ratio
copy_VarCoords(change,PR)


;-------------read model pd Variability-------------
variability := new((/Nvar,Nmodel,nlat,nlon/),float)

do ivar = 0,Nvar-1

  f1 := addfile(diri+"CMIP5/pr/"+var_name(ivar)+"/"+season+"/variability/multimodel/variability_"+var_name(ivar)+"_multimodel_"+season+"_1997-2014.nc","r")
  variability(ivar,0:N1-1,:,:) = f1->$VAR_measure$          ;(35model lat lon)

  f2 := addfile(diri+"CMIP6/pr/"+var_name(ivar)+"/"+season+"/variability/multimodel/variability_"+var_name(ivar)+"_multimodel_"+season+"_1997-2014.nc","r")
  variability(ivar,N1:N1+N2-1,:,:) = f2->$VAR_measure$       ;(25model lat lon)

end do


;----------------read model warming periods-----------------
warming_periods := new(Nmodel,integer)

f1 := addfile(diri+"CMIP5/tas/GMST/multimodel/GMST_multimodel_warming_levels_year.nc","r")
warming_periods(0:N1-1) = f1->year_center(:,ilevel)                   ;(35model)

f2 := addfile(diri+"CMIP6/tas/GMST/multimodel/GMST_multimodel_warming_levels_year.nc","r")
warming_periods(N1:N1+N2-1) = f2->year_center(:,ilevel)               ;(25model)


imodel_valid := ind(.not.(ismissing(warming_periods)))
Nmodel_valid := dimsizes(imodel_valid)


;**************read OBS variability******************
VAR_obs := new((/Nvar,Nobs,nlat,nlon/),float)

do ivar = 0,Nvar-1
  fobs := addfile(diri+"OBS/MultiOBS/"+var_name(ivar)+"/variability/variability_"+var_name(ivar)+"_MultiOBS_"+season+"_1997-2014.nc","r")
  VAR_obs(ivar,:,:,:) = fobs->$VAR_measure$               ;(4dataset lat lon)
end do



;***************************************
;mask to land
;***************************************
a := addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
lsdata := a->LSMASK

lsm := landsea_mask(lsdata,PR&latitude,PR&longitude)
lsm1 := conform(PR,lsm,(/3,4/))                 ;(/Nvar,Nmodel,Nindex,lat,lon/)

PR_land := where(((lsm1 .eq. 1) .or. (lsm1 .eq. 3)), PR, PR@_FillValue)
copy_VarCoords(PR,PR_land)

variability_land := where(((lsm1(:,:,0,:,:) .eq. 1) .or. (lsm1(:,:,0,:,:) .eq. 3)), variability, variability@_FillValue)
copy_VarCoords(variability,variability_land)

VAR_obs_land := where(((lsm1(:,0:Nobs-1,0,:,:) .eq. 1) .or. (lsm1(:,0:Nobs-1,0,:,:) .eq. 3)), VAR_obs, VAR_obs@_FillValue)
copy_VarCoords(VAR_obs,VAR_obs_land)


;***********************************************
;regional mean
;***********************************************
lat := variability&latitude
rad = 4*atan(1.0)/180.
clat := cos(lat*rad)
clat!0 = "lat"
clat&lat = lat

variability_regional := new((/Nvar,Nmodel,Nregion/),float)
PR_regional := new((/Nvar,Nmodel,Nindex,Nregion/),float)
VAR_obs_regional := new((/Nvar,Nobs,Nregion/),float)

do iregion = 0,Nregion-1
  variability_regional(:,:,iregion) = wgt_areaave_Wrap(variability(:,:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
  PR_regional(:,:,:,iregion) = wgt_areaave_Wrap(PR(:,:,:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
  VAR_obs_regional(:,:,iregion) = wgt_areaave_Wrap(VAR_obs(:,:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)

  if (domain(iregion).eq. "ocean") then
    VAR_obs_regional(:,1:2,iregion) = VAR_obs_regional@_FillValue
  end if            ;set GPCC & CPC to missing for ocean domains

end do

;-------------for land regions 4-5 ----------------
;subregions := (/"NH_extratropics","SH_extratropics","Northern_Asia","WNP_monsoon","Europe","Eastern_NAmerica"/)
;--------------------------------------------------
do iregion = 4,5
  variability_regional(:,:,iregion) = wgt_areaave_Wrap(variability_land(:,:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
  PR_regional(:,:,:,iregion) = wgt_areaave_Wrap(PR_land(:,:,:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
  VAR_obs_regional(:,:,iregion) = wgt_areaave_Wrap(VAR_obs_land(:,:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
end do




;-------------OBS uncertainty------------------

VAR_obs_regional_MME := dim_avg_n_Wrap(VAR_obs_regional,1)             ;(/Nvar,Nregion/)
;VAR_obs_regional_stddev := dim_stddev_n_Wrap(VAR_obs_regional,1)
;VAR_obs_regional_variance := dim_variance_n_Wrap(VAR_obs_regional,1)


;----------------------------------

corr := new((/Nvar,Nindex,Nregion/),float)               ;Pearson correlation coefficient
slope := new((/Nvar,Nindex,Nregion/),float)

PR_regional_corrected := PR_regional
PR_regional_corrected = PR_regional_corrected@_FillValue


do ivar = 0,Nvar-1
  do iindex = 0,Nindex-1
    do iregion = 0,Nregion-1

      slope(ivar,iindex,iregion) = regline( variability_regional(ivar,imodel_valid,iregion), PR_regional(ivar,imodel_valid,iindex,iregion) )
      corr(ivar,iindex,iregion) = escorc( variability_regional(ivar,imodel_valid,iregion), PR_regional(ivar,imodel_valid,iindex,iregion) )


      ;-------------remove linear bias------------
      do imodel = 0,Nmodel-1
        if (.not.(ismissing(warming_periods(imodel)))) then
          Xbias := variability_regional(ivar,imodel,iregion) - VAR_obs_regional_MME(ivar,iregion)
          Ybias := slope(ivar,iindex,iregion)*Xbias
          PR_regional_corrected(ivar,imodel,iindex,iregion) = PR_regional(ivar,imodel,iindex,iregion) - Ybias
        end if
      end do              ;{correct for each model}

    end do
  end do
end do


;============uncertainty range before & after correction=============
Yunc := PR_regional(:,imodel_valid,:,:)                 ;(/Nvar,model,Nindex,Nregion/)
Ycor := PR_regional_corrected(:,imodel_valid,:,:)

ii25 := toint(Nmodel_valid*0.10)
ii50 := toint(Nmodel_valid*0.50)
ii75 := toint(Nmodel_valid*0.90)

ip1 := dim_pqsort_n(Yunc, 2, 1)            ;increase & sort X
ip2 := dim_pqsort_n(Ycor, 2, 1)

Yunc_ii25 := Yunc(:,ii25,:,:)              ;(/Nvar,Nindex,Nregion/)
Yunc_ii50 := Yunc(:,ii50,:,:)
Yunc_ii75 := Yunc(:,ii75,:,:)

Ycor_ii25 := Ycor(:,ii25,:,:)
Ycor_ii50 := Ycor(:,ii50,:,:)
Ycor_ii75 := Ycor(:,ii75,:,:)

;-------------------print median estimates-------------------
print("median estimates: uncorrected      corrected")
iindex = 1    ;R95
ivar = 1      ;pr5d
do iregion = 0,3
  print( subregions(iregion)+": "+Yunc_ii50(ivar,iindex,iregion)+"    "+Ycor_ii50(ivar,iindex,iregion))
end do


;------------intermodel variance reduced----------------

variance_unc := dim_variance_n_Wrap(PR_regional(:,imodel_valid,:,:),1)               ;;(/Nvar,Nindex,Nregion/)
variance_cor := dim_variance_n_Wrap(PR_regional_corrected(:,imodel_valid,:,:),1)

variance_reduce := (variance_unc - variance_cor)/variance_unc*100.         ;(/Nvar,Nindex,Nregion/)
copy_VarCoords(PR_regional(:,0,:,:),variance_reduce)


;print(SNR(:,2)+"   "+corr(:,1,2))

;--------------PDF of PR------------------------
Nbin := 36
bin := fspan(0.5,4.0,Nbin)

bin(0) = 0.
bin(Nbin-1) = 10.
dbin = bin(5)-bin(4)

PR_pdf_unc := new((/Nvar,Nindex,Nregion,Nbin/),float)
PR_pdf_unc = 0.
PR_pdf_cor := PR_pdf_unc

do ivar = 0,Nvar-1
  do iindex = 0,Nindex-1
    do iregion = 0,Nregion-1
      do ibin = 0,Nbin-2
        ;-----------------uncorrected--------------
        idx1 := ind( (PR_regional(ivar,imodel_valid,iindex,iregion) .ge. bin(ibin)).and.(PR_regional(ivar,imodel_valid,iindex,iregion) .lt. bin(ibin+1)) )
        if .not.(all(ismissing(idx1))) then
          PR_pdf_unc(ivar,iindex,iregion,ibin) = dimsizes(idx1)
        end if

        ;------------------corrected---------------
        idx2 := ind( (PR_regional_corrected(ivar,imodel_valid,iindex,iregion) .ge. bin(ibin)).and.(PR_regional_corrected(ivar,imodel_valid,iindex,iregion) .lt. bin(ibin+1)) )
        if .not.(all(ismissing(idx2))) then
          PR_pdf_cor(ivar,iindex,iregion,ibin) = dimsizes(idx2)
        end if

      end do
    end do
  end do
end do

PR_pdf_unc = PR_pdf_unc/Nmodel_valid            ;occurrence frequency
PR_pdf_cor = PR_pdf_cor/Nmodel_valid




;***************************************************************
do iregion = 0,3           ;JJA 0-3 regions
do iindex = 1,1                      ;R95

wks := gsn_open_wks("pdf","/WORK/zhangwx/EConstraint/pic/correlation/onlyCMIP/correction/PDF_correction_slope_unadjusted_PR_regions_"+warming_levels(ilevel)+"_"+VAR_measure+"_"+index_name(iindex)+"_"+subregions(iregion)+"_"+season)

plot := new((/Nvar/),graphic)
plot_pl := plot
plot_xy_ref := plot
plot_pm := new((/Nvar,Nmodel/),graphic)


plot_tx := new((/Nvar,3/),graphic)

plot_unc := new((/Nvar/),graphic)
plot_cor := new((/Nvar/),graphic)
plot_unc_MME := new((/Nvar/),graphic)
plot_cor_MME := new((/Nvar/),graphic)

res = True
res@gsnFrame = False
res@gsnDraw = False

res@tmXTOn = False
res@tmYROn = False

res@gsnXYBarChart = True
res@gsnXYBarChartFillOpacityF = 0.3
res@gsnXYBarChartBarWidth = dbin/3.

;res@gsnXYBarChartOutlineOnly = True

;res@tmXBMajorOutwardLengthF = 0.
;res@tmXBMinorOutwardLengthF = 0.
;res@tmYLMajorOutwardLengthF = 0.
;res@tmYLMinorOutwardLengthF = 0.

res@tmBorderThicknessF = 1.0
res@tmXBMajorThicknessF = 1.0
res@tmXBMinorThicknessF = 1.0
res@tmYLMajorThicknessF = 1.0
res@tmYLMinorThicknessF = 1.0


;res@xyMarkLineModes = "Markers"
;res@xyMarkers =  16
;res@xyMarkerColor = "white"
;res@xyMarkerSizeF = 0.01

;res@tiMainFontHeightF = 0.034
res@tiXAxisFontHeightF = 0.023
res@tiYAxisFontHeightF = 0.023
res@tmXBLabelFontHeightF = 0.02
res@tmYLLabelFontHeightF = 0.02
res@gsnLeftStringFontHeightF = 0.028
res@gsnRightStringFontHeightF = 0.028
;res@gsnRightStringFontHeightF = 0.035
;res@gsnRightStringOrthogonalPosF = -0.15
;res@gsnRightStringParallelPosF = 0.93

res@gsnLeftString = ""
res@gsnRightString = ""

;res@xyLineThicknessF = 1.
;res@xyLineColor = "grey"

;-------------------------
res1 = res

colors := (/"grey","coral","orange"/)
colors_font := (/"grey30","sienna3","orange"/)

res@gsnXYBarChartColors = colors(0)
res@xyLineColors = colors(0)

res1@gsnXYBarChartColors = colors(1)
res1@xyLineColors = colors(1)


;-------------------------

plres = True
plres@gsLineThicknessF = 1.5
plres@gsLineColor = "grey70"
plres@gsLineDashPattern = 1



;-------------------------

markercolor = new(Nmodel,string)
markercolor(0:N1-1) = colors(0)
markercolor(N1:N1+N2-1) = colors(1)

modelcolor = new(Nmodel,string)
modelcolor(0:N1-1) = colors(0)
modelcolor(N1:N1+N2-1) = colors(1)

pmres = True
pmres@gsMarkerIndex =  16
pmres@gsMarkerSizeF = 0.005

;------------------------
txres = True
txres@txFontHeightF = 0.022
txres@txFontColor = "black"

txres_model = True
txres_model@txFontHeightF = 0.017


;*****************************************

leftst = (/"a","b","c","d","e","f"/)

;tmMin = new((/Nvar/),float)
;tmMax = tmMin
;height = tmMin

delta := dbin/6.

do ivar = 0,Nvar-1

  res@gsnLeftString = subregions_plot(iregion)
  res@gsnRightString = season

  res@tiXAxisString = "Probability ratio"
  res@tiYAxisString = "Frequency"

  res@trXMinF = min( (/ min(PR_regional(ivar,imodel_valid,iindex,iregion)), min(PR_regional_corrected(ivar,imodel_valid,iindex,iregion)) /))-0.1
  res@trXMaxF = max( (/ max(PR_regional(ivar,imodel_valid,iindex,iregion)), max(PR_regional_corrected(ivar,imodel_valid,iindex,iregion)) /))+0.1
  res@trYMinF = 0
  res@trYMaxF = max( (/ max(PR_pdf_unc(ivar,iindex,iregion,:)), max(PR_pdf_cor(ivar,iindex,iregion,:)) /) )*1.06

  plot_unc(ivar) = gsn_csm_xy (wks,bin-delta,PR_pdf_unc(ivar,iindex,iregion,:),res)              ;(/Nvar,Nindex,Nregion,Nbin/)
  plot_cor(ivar) = gsn_csm_xy (wks,bin+delta,PR_pdf_cor(ivar,iindex,iregion,:),res1)
  overlay(plot_unc(ivar),plot_cor(ivar))


  ;------------------MME------------------
  plres@gsLineColor = colors_font(0)
  plot_unc_MME(ivar) = gsn_add_polyline(wks,plot_unc(ivar),(/Yunc_ii50(ivar,iindex,iregion),Yunc_ii50(ivar,iindex,iregion)/),(/0,1/),plres)

  plres@gsLineColor = colors_font(1)
  plot_cor_MME(ivar) = gsn_add_polyline(wks,plot_unc(ivar),(/Ycor_ii50(ivar,iindex,iregion),Ycor_ii50(ivar,iindex,iregion)/),(/0,1/),plres)

  ;----------------------------------------
  ;refl := fspan(tmMin(ivar),tmMax(ivar),10)
  ;plot_xy_ref(ivar) = gsn_add_polyline(wks,plot(ivar),refl,refl,plres)

  ;----------------------------------------
  tx1 := sprintf("%4.2f",variance_reduce(ivar,iindex,iregion))+"%"
  txres@txFontColor = "black"
  plot_tx(ivar,2) = gsn_add_text(wks,plot_unc(ivar),tx1,res@trXMaxF-0.18,res@trYMaxF*0.9,txres)

  if (iregion .eq. 0) then
  txres@txFontColor = colors_font(0)
  plot_tx(ivar,0) = gsn_add_text(wks,plot_unc(ivar),"Uncorrected",res@trXMaxF-0.23,res@trYMaxF*0.82,txres)

  txres@txFontColor = colors_font(1)
  plot_tx(ivar,1) = gsn_add_text(wks,plot_unc(ivar),"Corrected",res@trXMaxF-0.2,res@trYMaxF*0.74,txres)
  end if


  ;-------------uncertainty range---------------
  ;-------------uncorrected-------------
  ;delta := height(ivar)*0.02
  ;plot_unc(ivar,0) = gsn_add_polyline(wks,plot(ivar),(/Yunc_ii25(ivar,iindex,iregion),Yunc_ii75(ivar,iindex,iregion)/),(/Ycor_ii50(ivar,iindex,iregion),Ycor_ii50(ivar,iindex,iregion)/), plres_un)
  ;plot_unc(ivar,1) = gsn_add_polyline(wks,plot(ivar),(/Yunc_ii25(ivar,iindex,iregion),Yunc_ii25(ivar,iindex,iregion)/),(/Ycor_ii50(ivar,iindex,iregion)-delta,Ycor_ii50(ivar,iindex,iregion)+delta/), plres_un)
  ;plot_unc(ivar,2) = gsn_add_polyline(wks,plot(ivar),(/Yunc_ii75(ivar,iindex,iregion),Yunc_ii75(ivar,iindex,iregion)/),(/Ycor_ii50(ivar,iindex,iregion)-delta,Ycor_ii50(ivar,iindex,iregion)+delta/), plres_un)

  ;------------corrected--------------
  ;plot_cor(ivar,0) = gsn_add_polyline(wks,plot(ivar),(/Yunc_ii50(ivar,iindex,iregion),Yunc_ii50(ivar,iindex,iregion)/),(/Ycor_ii25(ivar,iindex,iregion),Ycor_ii75(ivar,iindex,iregion)/), plres_un)
  ;plot_cor(ivar,1) = gsn_add_polyline(wks,plot(ivar),(/Yunc_ii50(ivar,iindex,iregion)-delta,Yunc_ii50(ivar,iindex,iregion)+delta/),(/Ycor_ii25(ivar,iindex,iregion),Ycor_ii25(ivar,iindex,iregion)/), plres_un)
  ;plot_cor(ivar,2) = gsn_add_polyline(wks,plot(ivar),(/Yunc_ii50(ivar,iindex,iregion)-delta,Yunc_ii50(ivar,iindex,iregion)+delta/),(/Ycor_ii75(ivar,iindex,iregion),Ycor_ii75(ivar,iindex,iregion)/), plres_un)



end do               ;{ivar}


;------------
txres2 = True
txres2@txFont = "helvetica-bold"
txres2@txFontHeightF = 0.012
gsn_text_ndc(wks,leftst(iregion),0.37,0.632,txres2)


;--------------------
resP = True
resP@gsnPanelTop = 0.95
;resP@gsnPanelLeft = 0.05
;resP@gsnPanelLabelBar  = True
;resP@gsnPanelYWhiteSpacePercent = 3.
resP@gsnPanelXWhiteSpacePercent = 3.
;resP@lbLabelFontHeightF = 0.012

gsn_panel(wks,plot_unc,(/1,3/),resP)


end do                  ;{iindex}
end do                  ;{iregion}


end
